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ABSTRACT 

The interstellar dust content in galaxies can be traced in extinction at optical wavelengths, or in emission in the far-infrared. Several 
studies have found that radiative transfer models that successfully explain the optical extinction in edge-on spiral galaxies generally 
underestimate the observed FIR/submm fluxes by a factor of about three. In order to investigate this so-called dust energy balance 
problem, we use two Milky Way-like galaxies produced by high-resolution hydrodynamical simulations. We create mock optical 
edge-on views of these simulated galaxies (using the radiative transfer code SKIRT), and we then fit the parameters of a basic spiral 
galaxy model to these images (using the fitting code FitSKIRT). The basic model includes smooth axisymmetric distributions along a 
Sersic bulge and exponential disc for the stars, and a second exponential disc for the dust. We find that the dust mass recovered by the 
fitted models is about three times smaller than the known dust mass of the hydrodynamical input models. This factor is in agreement 
with previous energy balance studies of real edge-on spiral galaxies. On the other hand, fitting the same basic model to less complex 
input models (e.g. a smooth exponential disc with a spiral perturbation or with random clumps), does recover the dust mass of the input 
model almost perfectly. Thus it seems that the complex asymmetries and the inhomogeneous structure of real and hydrodynamically 
simulated galaxies are a lot more efficient at hiding dust than the rather contrived geometries in typical quasi-analytical models. This 
effect may help explain the discrepancy between the dust emission predicted by radiative transfer models and the observed emission 
in energy balance studies for edge-on spiral galaxies. 
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1. Introduction 

Interstellar dust grains in galaxies play a special role in produc¬ 
ing and processing radiation. They efficiently absorb and scatter 
ultraviolet (UV), optical and near-infrared (NIR) photons and 
then reradiate the absorbed energy in the infrared and submm 
wavelength range. As a result, the presence of dust influences 
the apparent view of a galaxy and apparent structural param¬ 
eters such as scale lengths, surface brightnesses, luminosities 
and colours. The effects of dust attenuation on these parameters 
are nontrivial and sometimes counterintuitive (Byun et al. 1994; 
Tuffs et al. 2004; Pierini et al. 2004; Gadotti et al. 2010; Pastrav 
et al. 2013a,b). A major source of complexity is the relative ge¬ 
ometry between stars and dust in a galaxy: the same amount of 
dust can give rise to completely diflTerent levels of attenuation 
depending on the star-dust geometry (Witt et al. 1992; Baes & 
Dejonghe 2001b). Before we can correct for the attenuating ef¬ 
fects in a galaxy, we need to understand the total amount and 
spatial distribution of interstellar dust. 

Broadly speaking, the dust in galaxies can be traced in two 
ways. The first way is by studying the thermal emission of the 
dust at far-infrared (FIR) and submm wavelengths. We can deter¬ 
mine the total dust mass of a galaxy with a reasonable accuracy 
by fitting a modified blackbody or more advanced models to the 
FIR/submm observed spectral energy distribution. Especially af¬ 
ter the launch of the Herschel Space Observatory (Pilbratt et al. 
2010), a wealth of FIR/submm data has become available that 
allowed different teams to determine dust masses in thousands 


of galaxies (e.g., Dunne et al. 2011; Dale et al. 2012; Smith 
et al. 2013; Ciesla et al. 2014). Due to the diffraction limit in 
the FIR/submm and the limited diameter of space observatories, 
detailed investigations of the distribution of the dust remains dif¬ 
ficult, except for the most nearby galaxies (e.g., Kramer et al. 
2010; Fritz et al. 2012). 

An alternative method, which does not suffer from the poor 
spatial resolution affecting FIR/submm observations, consists of 
carefully modelling the attenuation effects of the dust on the 
stellar emission in the optical window using advanced radia¬ 
tive transfer techniques. Today, a range of powerful 3D radia¬ 
tive transfer codes is available (for an overview, see Steinacker 
et al. 2013), and the fitting of the radiative transfer models to data 
is being optimised beyond the elementary chi-by-eye (Xilouris 
et al. 1997; Steinacker et al. 2005; Bianchi 2007; Schechtman- 
Rook et al. 2012; De Geyter et al. 2013, 2014). 

Traditionally, detailed radiative modelling has been applied 
mainly to edge-on spiral galaxies, as the dust is then clearly vis¬ 
ible as prominent dust lanes in optical images. Most studies fo¬ 
cus on a single galaxy (e.g., Kylafis & Bahcall 1987; Xilouris 
et al. 1997, 1998; Popescu et al. 2000, 2011; Baes et al. 2010; 
De Fooze et al. 2012a,b; Schechtman-Rook et al. 2012) or mod¬ 
est sets of edge-on spiral galaxies (Xilouris et al. 1999; Bianchi 
2007; MacFachlan et al. 2011; De Geyter et al. 2014). 

A comparison of the extinction and the FIR/submm emis¬ 
sion, i.e. a study of the dust energy balance, gives the strongest 
constraints on the dust content of spiral galaxies. Dust energy 
balance studies of edge-on spiral galaxies reveal a discrepancy 
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between the FIR/submm emission predicted by radiative transfer 
models and the observed emission. Although the radiative trans¬ 
fer models explain the optical extinction very well, they typi¬ 
cally underestimate the observed dust emission by a factor of 
about three (Popescu et al. 2000; Misiriotis et al. 2001; Alton 
et al. 2004; Dasyra et al. 2005; Baes et al. 2010; De Looze et al. 
2012a,b). This has become generally known as the dust energy 
balance problem in edge-on spiral galaxies. 

Two broad scenarios have been proposed to explain this 
problem. A first possibility is that FIR/submm emissivity, stan¬ 
dardly taken from semi-empirical dust models based on the in¬ 
terstellar medium (ISM) in the Milky Way, is underestimated 
by a factor of a few. This idea was advocated by Alton et al. 
(2004) and Dasyra et al. (2005), and is corroborated by the large 
range of empirical values for emissivity circulating in the liter¬ 
ature. However, based on a detailed study of the edge-on spiral 
UGC 4754, Baes et al. (2010) argued that this possibility cannot 
explain the discrepancy. Indeed, they found an incompatibility 
between the absorbed stellar luminosity as obtained from the ra¬ 
diative model fits to the optical data and the emitted dust lumi¬ 
nosity in the FIR/submm that cannot be lifted by just increasing 
the value of the dust emissivity. 

The alternative scenario to explain the dust energy balance 
problem in spiral galaxies is geometrical in nature. As already 
indicated, the relative star-dust geometry in a galaxy is extremely 
important for the amount of attenuation. If a sizeable fraction of 
the dust is distributed in such a way that it hardly attenuates the 
bulk of the starlight, one could easily underestimate the amount 
of dust from radiative transfer modelling. At the same time, this 
dust can still contribute substantially to the observed FIR/submm 
emission. How exactly this additional dust would have to be dis¬ 
tributed is subject to debate, and various options have been pro¬ 
posed. 

One option is a very thin and dense dust disc next to the 
thicker dust disc that is responsible for the dust lane in the optical 
(Popescu et al. 2000, 2011; Misiriotis et al. 2001; Driver et al. 
2007). For example, Schechtman-Rook & Bershady (2014) find 
a very thin dust disk in certain edge-on spiral galaxies. On the 
other hand, Dasyra et al. (2005) showed that such a disc would 
have an observational signature at NIR wavelengths, which is 
at odds with the observations of the prototypical edge-on spiral 
galaxy NGC891. 

A second option is that most of the dust is locked up in so- 
called “clumps”, a rather vaguely defined term used in the radia¬ 
tive transfer community to indicate any form of inhomogeneities 
relative to a smooth medium. One can consider large-scale in¬ 
homogeneities such as bars and spiral arms, as well as small- 
scale structures such as dusty molecular clouds with or without 
embedded young stars. Thanks to the advancement of 3D dust 
radiative transfer, the effect of a non-homogeneous multi-phase 
dusty medium has been investigated by various teams (e.g., Witt 
& Gordon 1996, 2000; Wolf et al. 1998; Varosi & Dwek 1999; 
Bianchi et al. 2000; Hegmann & Kegel 2003; Indebetouw et al. 
2006). While the dusty ISM in spiral galaxies is far from smooth 
and homogeneous, the radiative transfer models being fitted to 
the optical images in energy balance studies are usually smooth 
and axisymmetric. This simplification may very well affect the 
results of these studies. 

Whether or not an inhomogeneous distribution of the dust, 
on small and/or on large scales, is a possible solution for the 
dust energy balance problem is still an open question. Several 
studies have attempted to address this issue. Concerning the role 
of spiral arms, the work of Misiriotis et al. (2000, hereafter MOO) 
is the most advanced study. They set up a suite of galaxy models 


with an analytical spiral structure, created mock edge-on optical 
images for these galaxies, and subsequently fitted these simu¬ 
lated observations with a smooth, axisymmetric radiative trans¬ 
fer model. The fitted models were surprisingly accurate; essen¬ 
tially all parameters of the input model (dust mass, inclination, 
scale parameters of stars and dust, etc.) were properly recovered. 
In a similar way, Misiriotis & Bianchi (2002, hereafter M02) in¬ 
vestigated the effect of clumps or small-scale inhomogeneities: 
rather than models with a spiral perturbation they adopted the 
clumpy disk galaxy models of Bianchi et al. (2000) as input 
models. They concluded that, for the range of clumpy distribu¬ 
tions they considered, the neglect of clumping results in a sys¬ 
tematic underestimate of the dust mass. The underestimate was 
found to be never larger than 40%, however. 

These results seem to cast doubts on the ability of inhomo¬ 
geneities such as spiral arms or large molecular clouds to pro¬ 
vide an answer to the dust energy balance problem. However, 
one needs to take into account that the input models used in 
the studies by MOO and M02 were still relatively well-behaved. 
They were constructed using theoretical perturbations on essen¬ 
tially the same exponential discs as those used in the fitting, and 
still featured a rather symmetric and regular geometry. They are 
therefore a relatively poor representation for real galaxies, which 
are generally characterised by a much larger degree of geomet¬ 
rical complexity and asymmetry. 

The goal of this paper is to investigate whether a complex 
and inhomogeneous dusty ISM could provide an answer to the 
dust energy balance problem. In Sect. 2 we start with basic input 
models similar to those used by MOO and M02, and we verify 
that our fitting procedure recovers similar results. In Sect. 3 we 
apply the same approach to galaxy snapshots produced by state- 
of-the-art hydrodynamic simulations. We consider two models 
from different simulations, both of which show a complex and 
asymmetric geometry similar to our own Milky Way. In Sect. 4 
we discuss and interpret the results, and a summary is presented 
in Sect. 5. 

2. Basic input modeis 

Before we set out to model mock images from hydrodynamical 
simulations and compare the results to previous work, we first 
need to investigate whether we can reproduce the results of this 
previous work. In particular, we should verify that our fitting 
routine can reproduce the results of MOO and M02 for basic in¬ 
put models. This is not as obvious as it may seem. Indeed, while 
we follow a similar approach, there are potentially relevant dif¬ 
ferences in the modelling mechanisms. For example, the radia¬ 
tive transfer code used by MOO and M02 adopts the so-called 
method of scattered intensities (for details, see Kylafis & Bahcall 
1987; Xilouris et al. 1997; Baes & Dejonghe 2001a), whereas the 
fundamental algorithm in our SKIRT code is the Monte Carlo 
method. Perhaps more importantly, while MOO and M02 fit each 
model to a single image, we use oligochromatic fitting, i.e. we fit 
a single model simultaneously to images in the five SDSS optical 
bands. This approach can help to eliminate degeneracies in the 
parameter space (De Geyter et al. 2014). Hence, we set up two 
test cases that are very similar to the basic input models explored 
by MOO and M02. 

2.1. The input models 

The first step in our modelling is the setup of the models, i.e., 
the choice of the basic model from which to start and the per¬ 
turbations that we apply this model. For the underlying model. 
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we assume a smooth and axisymmetric model, consisting of a 
double-exponential stellar disc, a flattened Sersic bulge, and a 
double-exponential dust disc (see also, e.g., Xilouris et al. 1999; 
Bianchi 2007; MacLachlan et al. 2011). For the parameters of 
the underlying model, we use the average values obtained by fit¬ 
ting 10 real galaxies with FitSKIRT from the CALIFA survey 
(see Table 4 in De Geyter et al. 2014). 

First, we consider a spiral structure perturbation, similar to 
the models by MOO. Both the stars in the disk and the dust are 
perturbed into a logarithmic spiral arm pattern. For the parame¬ 
ters of the spiral arm perturbation, we select an average model 
also considered in MOO. There are two spiral arms, the weight of 
the spiral perturbation is 30% for the stars and 40% for the dust, 
and the spiral arm pitch angle is 20 degrees. 

Secondly, we set up an idealised model for a clumpy disc 
galaxy similar to M02. In this case, we divide the dusty medium 
into a smooth and a clumpy component. We take a difiTerent ap¬ 
proach from Bianchi et al. (2000), who adopted the two-phase 
medium algorithm explored by many authors (Witt & Gordon 
1996, 2000; Matthews & Wood 2001; Wolf et al. 1998; Stalevski 
et al. 2012). Instead, we deposit half of the total dust mass in 
10,000 individual spherical clumps, each with an outer radius of 
300 pc and an internal cubic spline kernel distribution (Hockney 
& Eastwood 1981). Each of these clumps is positioned at a ran¬ 
dom location chosen according to the underlying smooth dust 
density field (see also De Looze et al. 2014; Camps & Baes 
2015). 

2.2. Creation of mock images 

The second step in our analysis is the creation of mock images 
for these two models. In particular, we need images in the SDSS 
ugriz bands and in the edge-on orientation, so that they can be 
used subsequently for radiative transfer modelling. The mock 
images were created with SKIRT, a 3D Monte Carlo dust ra¬ 
diative transfer code designed to simulate continuum radiation 
transfer in dusty astrophysical systems (Baes et al. 2003, 2011; 
Camps & Baes 2015). It supports multiple anisotropic scatter¬ 
ing, absorption and re-emission by interstellar dust (including 
stochastically heated grains), and can create spectral energy dis¬ 
tributions and images at arbitrary wavelengths and from arbitrary 
points of view. SKIRT is publicly available from a GitHub code 
repository.^ 

2.3. Radiative transfer modelling 

Einally, the third step in the modelling consists of fitting radiative 
transfer models to the mock images, as is done for real galaxies. 
We use the EitSKIRT code for this goal. EitSKIRT (De Geyter 
et al. 2013, 2014) is a tool designed to recover the 3D spatial 
distribution of stars and dust by fitting radiative transfer models 
to optical images. The code reads any number of images in dif¬ 
ferent bands and simultaneously fits a radiative transfer model to 
all images. The model can include an arbitrary combination of 
stellar and dust components. It combines SKIRT with the power 
of the genetic algorithms-based optimisation library GAlib (Wall 
1996) to perform the actual fitting. 

In this work, we use the same radiative transfer fitting model 
as the one employed by De Geyter et al. (2014) to fit real edge- 
on spiral galaxies. We fit a smooth axisymmetric model, similar 
to the underlying model, to the mock images of the perturbed 

^ SKIRT code repository: https://github.coin/skirt/skirt, 
SKIRT documentation: http: //www. skirt. ugent. be 


Table 1. Model parameter values recovered by the EitSKIRT ra¬ 
diative transfer fits for two basic input models. The third column 
provides the “true” values of the input models, and the fourth 
and fifth columns list the values recovered for the model with a 
spiral arm perturbation and a clumpy structure respectively. Eor 
a definition of each parameter, see De Geyter et al. (2014). 


Parameter 

Units 

Input 

Spiral 

Clumpy 

Hr,. 

kpc 

4.23 

4.26 + 0.03 

4.48 + 0.06 

K,* 

kpc 

0.51 

0.59 + 0.03 

0.58 + 0.01 

^efF 

kpc 

2.31 

2.44 + 0.15 

2.52 + 0.09 

n 

— 

2.61 

2.27 + 0.12 

2.0+ 0.1 

q 

— 

0.56 

0.58 + 0.01 

0.55 + 0.01 


kpc 

6.03 

4.89 + 0.56 

4.94 + 0.22 

Kd 

kpc 

0.23 

0.201 + 0.004 

0.203 + 0.003 

Md 

10^ Mo 

3.02 

2.34 + 0.12 

2.5+ 0.1 

i 

deg 

90 

90.01+0.01 

89.96 + 0.02 


model. The radiative transfer modelling essentially comes down 
to a strongly nonlinear ;^^-minimalisation in a 21-dimensional 
parameter space (we consider five parameters describing the 
stellar geometry, two parameters for the dust geometry, the lu¬ 
minosity of the bulge and disc in each of the five bands, the total 
dust mass and three projection parameters). Eor more details, we 
refer to Section 3.1 of De Geyter et al. (2014). 

The fitting results are presented in Eig. 1 and Table 1. The 
leftmost column in Eig. 1 represents the original mock images in 
the five bands, the central column shows the images correspond¬ 
ing to the best fitting model, and the rightmost column presents 
the residual images, expressed as the relative difference between 
the surface brightness of the real image and the model fit. The 
residual frames for the first input model show some discrepan¬ 
cies to the left of the centre of the galaxy, corresponding to a 
maximum in the spiral arm perturbation. Even in this feature, 
the relative difference between model and fit is only of the order 
of 20%, so this is an excellent fit. Eor the clumpy input model, 
the largest discrepancies occur in the u and g bands, where the 
effects of the dust are most prominent, but in general the quality 
of the fit is very satisfactory as well. 

Table 1 compares the parameters recovered by the fitting pro¬ 
cedure to the corresponding parameters of the input models. Eor 
both models, all input parameters are recovered to within 25%, 
and often even better. Specifically, the total dust mass is under¬ 
estimated by less than 25% in both cases. These results are in 
line with those reported by MOO and M02. 

3. Input models from hydrodynamical simulations 

Now that we have shown that we can reproduce the previous 
results by MOO and M02, we can move to the next level, and 
consider more realistic models for spiral galaxies. The modelling 
follows the same strategy as we used for the basic input models 
considered in the previous section. 

3.1. The input models 

In recent years, hydrodynamical simulations have started to suc¬ 
cessfully reproduce late-type spiral galaxies (Boumaud et al. 
2007, 2009; Govemato et al. 2009; Agertz et al. 2011; Guedes 
et al. 2011; Wada et al. 2011; Stinson et al. 2013; Renaud et al. 
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Fig. 1. Results of the FitSKIRT radiative transfer fits for two basic input models; one with analytical spiral perturbation (upper 
half) and one with random clumps (lower half). The left column shows the reference images produced by SKIRT in each of the u, 
g, r, /, and z bands. The middle column shows the corresponding fit obtained with FitSKIRT. The right column contains residual 
images showing the relative deviation between the fit and the reference image. The colour bar at the bottom presents the scale of the 
deviation in the residual images. 


2013; Marinacci et al. 2014; Inoue & Saitoh 2014). The spatial 
resolution of these simulations is sufficient to resolve both large- 
and small-scale inhomogeneities. We consider snapshots from 
two different simulations for the analysis in this paper. 

The first input model is taken from the simulation by Renaud 
et al. (2013), hereafter called the R13 simulation. This is a 
self-consistent hydrodynamical simulation of a Milky Way-like 
galaxy performed with the Adaptive Mesh Refinement (AMR) 
code RAMSES (Teyssier 2002). The goal of this simulation was 
to reproduce an isolated grand-design spiral galaxy and to fo¬ 


cus on the structure of the interstellar medium at the highest 
resolution possible. The simulation resolves the ISM down to 
scales of 0.05 pc, and includes stellar feedback through photo¬ 
ionization, radiative pressure and supernovae. The simulation 
started with 30 million particles in the dark matter halo and an¬ 
other 30 million particles representing the stars distributed over 
a bulge, spheroid, thin disc and thick disc. The gaseous disc is 
represented with 240 million AMR cells varying in size depend¬ 
ing on the structure of the ISM. The mass resolution is 160 M© 
for young star particles and 30 M© for the gas in the most refined 
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Fig. 2. A cut through the dust density distribution along the equatorial plane of the R13 galaxy (left panel) and of the Eris galaxy 
(right panel). The dust density is given in M© pc“^. 


cells. The initial components in the simulation are axisymmetric, 
but non-axisymmetric structures such as a prominent bar and spi¬ 
ral structure are developed quickly because of the instabilities in 
the velocity distribution. 

We use the snapshot corresponding to the final state of the 
galaxy, after a run time of 800 Myr. We consider a box with di¬ 
mensions 20 X 20 X 4 kpc^, which contains a total stellar mass of 
4.4 X 10^^ Mq and a total gas mass of 3.0 x 10^ M©. The galaxy 
hosts a prominent bulge and spiral arms, and is characterised by 
a high star formation rate (SFR); the total SFR for stars younger 
than 57 Myr is about 7.3 M© yr“^. The bar has a central area 
of 1 kpc^ that contains almost no ongoing star formation (see 
Emsellem et al. 2015), whereas dense star-forming clouds are 
abundant at the outer regions of the bar. Along the spiral arms, 
dense clumps of gas and star forming regions have formed with a 
relatively uniform separation (the so-called “beads on a string”). 
This high level of concentration is related to the relatively short 
simulation run time; continuing the simulation would cause stel¬ 
lar feedback to more evenly spread matter along the disc. 

The second input model in our study is taken from the Eris 
simulation (Guedes et al. 2011), one of the most advanced and 
realistic simulations of the formation of a Milky Way class 
galaxy. Eris was set up as a zoom-in cosmological simulation, 
powered by the N-body/SPH GASOLINE code (Wadsley et al. 
2004). The simulation follows the formation of a galaxy halo 
with mass Myir = 7.9x 10^^ M© from z = 90 to the present epoch 
in a full cosmological setting. The target halo is sampled with 26 
millions particles divided equally between the dark matter par¬ 
ticles and gas particles. Apart from the obvious gravity and hy- 
drodynamical forces, the simulation includes Compton cooling, 
atomic cooling, metallicity-dependent radiative cooling at low 
temperatures, the ionising effect of a uniform UV background, 
star formation and supernova feedback. 


We use the final snapshot corresponding to the present 
epoch. At z = 0, Eris is a Milky Way-like galaxy characterised 
by a beautiful spiral structure and a small bulge in the centre. 
The structural properties, the mass budget in the various com¬ 
ponents, and the scaling relations between mass and luminosity 
are consistent with a host of observational constraints. For our 
analysis, we consider all particles in a box of 28 x 28 x 6 kpc^. 
This box has 7.8 million star particles and 0.25 million gas par¬ 
ticles, with a total stellar mass of 3.5 x 10^^ M© and a total gas 
mass of 4.3 X 10^ M© (about 5 x 10^ M© per star particle and 
2 X 10"^ Mq per gas particle). The smoothing length for the gas 
particles ranges between 56 and 2455 pc. 

3.2. Creation of mock images 

The second step in the analysis is again the creation of mock 
images. In order to do so, SKIRT was extended with modules 
that read the output of hydrodynamical simulations as input for 
the radiative transfer calculations (see also Schaye et al. 2015). 
We used the GALAXEV library of simple stellar populations 
(Bruzual & Chariot 2003) to determine the intrinsic emission of 
the stars in the simulation. For the optical properties of the dust 
we use the BARE-GR-S dust model from Zubko et al. (2004), 
which is fine-tuned to be consistent with the extinction, diffuse 
emission and depletion in the Milky Way. 

One particular aspect that needs special care is the determi¬ 
nation of the dust density. As the hydrodynamical simulations do 
not track the dust directly, a recipe needs to be chosen to set the 
dust density from the properties of the gaseous medium. We use 
the assumption that a constant fraction of the metals in the ISM 
is locked up into dust grains. While observations show metallic- 
ity gradients across galaxies (Dobashi et al. 2008; Paradis et al. 
2012) and local variations in the fraction of metals locked into 
dust grains (Hirashita 2012), this simplification does not affect 
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Fig. 3. Mock face-on views for the R13 galaxy (left panel) and for the Eris galaxy (right panel). The three-colour images are based 
on the r, g and u band images produced by the SKIRT radiative transfer code. 



Fig. 4. Mock edge-on views for the R13 galaxy, looking at the head of the bar (left panel), and for the Eris galaxy (right panel). The 
three-colour images are based on the r, g and u band images produced by the SKIRT radiative transfer code. 


our analysis, because we only need to create mock images re¬ 
flecting a certain dust mass. In particular, we set the 3D dust 
density using 

Pdust — fdust ^ Pgas (1) 

where Z is the metallicity of the gas, and /dust is the fraction of 
metals locked up in dust grains. Eor the R13 simulation we use 
/dust = 0-3, the same value as used in the EAGLE simulation 
(Schaye et al. 2015). Eor the Eris simulation, we use a slightly 
larger value, /dust = 0-5, since this galaxy has a relatively low 
metal content. Both values are within the range suggested by 
observations, roughly between 0.2 and 0.7 (Dwek 1998; James 
et al. 2002; Mattsson & Andersen 2012; De Cia et al. 2013; Zafar 
& Watson 2013). Eor both simulations, the gas density can be 
determined from the simulation snapshot. In the R13 simulation, 
the metallicity is not stored or evolved, and we adopt a flxed solar 
metallicity. Combined with the flxed value for /dust this comes 
down to a constant gas-to-dust ratio of 166, which is in good 
agreement with the gas-to-dust ratio in the Milky Way and other 
spiral galaxies (Zubko et al. 2004; Draine et al. 2007; Remy- 
Ruyer et al. 2014). Eor the Eris simulation, the gas metallicity 
is stored for every SPH gas particle, such that the dust density 
can be calculated. Using a higher value of /dust boosts the total 
dust mass, even with the rather low metallicities in this galaxy 
(ranging from zero up to Z = 0.266). 


Another crucial ingredient for the present simulations is the 
setup of the grid on which the dust density is defined. SKIRT 
can handle any 3D geometry, thanks to the efficient use of hier¬ 
archical and unstructured grids for the dust medium (Saftly et al. 
2013, 2014; Camps et al. 2013). Eor the simulations in this pa¬ 
per, we used a hierarchical k-d tree with a resolution up to 0.6 pc 
for the R13 and 6.8 pc for the Eris simulation. Eor both simula¬ 
tions, the final grid contains about 1.5 million dust cells. Eigure 2 
shows the dust density in the equatorial plane for both models. 

Eor each of the two input models, we created mock images 
in the SDSS ugriz bands, in face-on and edge-on configurations 
(although we only need the edge-on views for our study). Three- 
colour images based on the r, g and u bands are shown in Eigs. 3 
and 4. In spite of the complex intrinsic structure of both the stel¬ 
lar and dust distribution apparent in the face-on views, the edge- 
on morphology seems rather smooth and regular. 

3.3. Radiative transfer modelling 

Einally, we fit radiative transfer models to the mock images using 
EitSKIRT, where we use exactly the same approach as for the 
basic input models from the previous section. The results for 
both snapshots are shown in Eig. 5 and Table 2. 
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Fig. 5. Results of the FitSKIRT radiative transfer fits for the R13 galaxy, looking at the head of the bar (upper half), and for the 
Eris galaxy (lower half). The left column shows the reference images produced by SKIRT in each of the w, g, r, /, and z bands. 
The middle column shows the corresponding fit obtained with FitSKIRT. The right column contains residual images showing the 
relative deviation between the fit and the reference image. The colour bar at the bottom presents the scale of the deviation in the 
residual images. 


The global morphology of each input model is accurately 
reproduced in the fit, although some features are not fitted com¬ 
pletely. In particular, the residual frames for the R13 model show 
some discrepancies in the dust lane area. This structure is most 
strongly present in the shortest wavelength bands, but it per¬ 
sists up to the z band. These discrepancies correspond to clumpy 
structures, both in the stellar distribution (the beads on a string) 
and in the dust distribution. However, most of the pixels in the 
residual frames have a discrepancy of less than 30%. In fact, the 
central areas of the residual frames for the R13 model shows 
an obvious similarity to those corresponding to the clumpy disc 


model (Fig. 1, lower panels). For the Eris simulation, the dust 
lane is not very prominent and it has a discontinuous shape 
with a lot of clumpy and irregular structures up to the edges 
of the galaxy, which makes it a hard galaxy to fit. The resid¬ 
ual frames show discrepancies around the bulge corresponding 
to the molecular clouds and star forming regions. 

Table 2 lists the most important model parameter values re¬ 
covered by the radiative transfer fits for each galaxy. To gauge 
the effect of the large-scale asymmetry in the R13 galaxy due 
to its prominent bar, we performed two independent fitting pro¬ 
cedures based on perpendicular edge-on viewpoints, respec- 


7 





















W. Saftly et al.: Dust energy balance problem 


Table 2. Model parameter values recovered by the FitSKIRT ra¬ 
diative transfer fits for the R13 and Eris galaxies. For the R13 
galaxy, we list two sets of parameters (head and side). They cor¬ 
respond to two independent radiative transfer fits, based on per¬ 
pendicular edge-on viewpoints, respectively looking at the head 
and the side of the bar. For a definition of each parameter, see 
De Geyter et al. (2014). 


Par. Units R13 (head) R13 (side) Eris 


hR,* 

kpc 

4.17 

+ 

0.38 

3.57 

+ 

0.39 

4.14 

+ 

0.13 

K, 

kpc 

0.43 

+ 

0.01 

0.40 

+ 

0.01 

0.37 

+ 

0.01 

^eff 

kpc 

0.93 

+ 

0.04 

1.12 

+ 

0.08 

0.49 

+ 

0.07 

n 

— 

0.64 

+ 

0.03 

0.80 

+ 

0.04 

6.71 

+ 

0.48 

q 

— 

0.47 

+ 

0.01 

0.40 

+ 

0.03 

0.56 

+ 

0.03 


kpc 

9.47 

+ 

0.69 

4.23 

+ 

0.80 

0.35 

+ 

0.05 

Kd 

kpc 

0.064 

+ 

0.005 

0.091 

+ 

0.006 

0.095 

+ 

0.008 

Md 

10*^ Me 

6.26 

+ 

0.24 

7.68 

+ 

0.33 

1.48 

+ 

0.18 

i 

deg 

89.96 

+ 

0.02 

89.96 

+ 

0.02 

90.04 

+ 

0.16 


tively looking at the head and the side of the bar. Given the 
dust masses of the input models (Md,Ri 3 = 1.81 x 10^ M©; 
^d,Eris = 5.94 X 10^ M©) and the recovered values in Table 2, 
it follows that the radiative transfer model underestimates the 
“true” dust mass by a factor of 2.9 + 0.2 respectively 2.4 + 0.2 
for the R13 galaxy, and by a factor of 4.0+0.4 for the Eris galaxy. 

4. Discussion 

The previous work by MOO and M02 cited in Sect. 1, and our 
results presented in Sect. 2, indicate that fitting a smooth galaxy 
model to a certain class of basic, quasi-analytical input models 
recovers the intrinsic dust mass of the input model to within 40% 
or better. In fact, as evidenced by the close correspondence be¬ 
tween recovered and input values listed in Table 1 , a logarithmic 
spiral arm perturbation or a clumpy dust distribution seem to 
have only a modest effect on the structural parameters observed 
in an optical edge-on view of a disc galaxy, including the derived 
dust mass. 

In contrast, our results presented in Sect. 3 indicate that fit¬ 
ting a smooth galaxy model to more realistic galaxy models, ob¬ 
tained from high-resolution hydrodynamical simulations, under¬ 
estimates the intrinsic dust mass of the input model by a factor of 
about three. This is a tantalising result, especially since a factor 
of roughly the same magnitude has been found in energy bal¬ 
ance studies of real edge-on spiral galaxies (Popescu et al. 2000; 
Misiriotis et al. 2001; Alton et al. 2004; Dasyra et al. 2005; Baes 
et al. 2010; De Looze et al. 2012a,b). 

It is tempting to conclude that the level of dust underesti¬ 
mation is driven by the fundamental differences between the in¬ 
put models. The models in MOO and M02, and in our Sect. 2, 
are derived from well-behaved, smooth disc models by apply¬ 
ing a relatively modest perturbation. For example, the spiral arm 
perturbation is fully analytical and cancels out exactly when av¬ 
eraged over azimuth. The galaxy models constructed from hy¬ 
drodynamical simulation snapshots, presented in Sect. 3, feature 
much more realistic inhomogeneities at a wide range of scales, 
from large-scale bars and spiral arms to parsec-sized clumps and 
filaments. These structural complexities may very well be re¬ 
sponsible for a higher level of dust underestimation in the radia¬ 
tive transfer fits. 


In other words, our modelling suggests that the complex and 
inhomogeneous structure of galaxies can hide up to three times 
more dust than is “observed” when the optical images are fitted 
with smooth axisymmetric models. FIR/submm observations of 
several spiral galaxies also imply a factor of three times more 
dust than visible in the optical; this correspondence suggest that 
the inhomogeneous structure of the ISM possibly is the source 
of the dust energy balance problem. The recent work by De 
Looze et al. (2014) supports this hypothesis. They performed 
a detailed panchromatic radiative transfer modelling of the face- 
on galaxy M51 with a model that includes the complex geom¬ 
etry as derived from the FUV attenuation map. The model self- 
consistently reproduces the surface brightness images from UV 
to submm wavelengths. The face-on analysis is of course less af¬ 
fected by optical depth effects along the line of sight, which may 
have contributed to this result as well. 

We must be careful not to jump to conclusions. First, while a 
typical factor of about three is found by several teams for differ¬ 
ent edge-on spiral galaxies (Popescu et al. 2000; Misiriotis et al. 
2001; Alton et al. 2004; Dasyra et al. 2005; Baes et al. 2010; 
De Looze et al. 2012a,b), this is by no means an ubiquitous fea¬ 
ture. This was shown most recently by De Geyter et al. (2015), 
who performed the same fitting procedure as presented in this 
paper on two edge-on spiral galaxies from the sample analysed 
in De Geyter et al. (2014). For one of the two galaxies, a typical 
factor-of-three discrepancy is observed between the best fitting 
FitSKIRT model and the observed FIR/submm SED, whereas 
for the other galaxy the FitSKIRT model accurately describes 
the observed spectrum both in absolute values and shape. 

Secondly, even for those galaxies in which a dust energy bal¬ 
ance problem is encountered, it is up to debate whether this can 
be ascribed to the same physical scenario. For example, the dust 
emission excess in the Sombrero galaxy is shown to be compati¬ 
ble with an additional unresolved cold dust reservoir (De Looze 
et al. 2012b), whereas the excess emission in the edge-on spiral 
UGC 4754 is rather compatible with an additional warmer com¬ 
ponent, such as expected when linked to recent star formation 
(Baes etal. 2010, 2011). 

Finally, our present study is based on just two simulated spi¬ 
ral galaxies, and each of them has its strengths and weaknesses. 
In both galaxies we recognise the structures and morphologies 
of real galaxies, including spiral arms, bars, bulges, star form¬ 
ing regions, and compact clumps. However, the R13 galaxy con¬ 
tains star forming regions that look somewhat artificial, and the 
galaxy’s dust lane is very thin and extended; and the Eris galaxy 
has a faint and rather fuzzy dust lane that is not visible in the r, i 
and z bands, which is atypical for real galaxies. 

5. Conclusion 

We set out to shed light on the dust energy balance problem in 
edge-on spiral galaxies by performing the radiative transfer fit¬ 
ting procedure that has been used previously for studying real 
galaxies on snapshots obtained from hydrodynamical simula¬ 
tions. These simulated galaxies feature a more realistic inhomo¬ 
geneous structure than typical quasi-analytical models, and their 
“true” dust mass is known. 

We used two simulated Milky Way-like galaxies as input 
models. The R13 simulation is a self-consistent hydrodynami¬ 
cal simulation performed with the AMR code RAMSES. The 
Eris simulation is a zoom-in cosmological simulation performed 
by the N-body/SPH GASOLINE code. We used our radiative 
transfer code SKIRT to create mock observational images in the 
SDSS ugriz bands for an edge-on view of both galaxies, and we 
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fitted the parameters of a basic, smooth disc galaxy model to 
these images with our radiative transfer fitting code FitSKIRT. 
We found that, for both galaxies, the dust mass is underesti¬ 
mated by a factor of about three. This result is strikingly close 
to what has been found in previous work by several teams who 
performed similar analyses for real edge-on spiral galaxies. 

In contrast, previous studies have shown that fitting a smooth 
disc galaxy model to a modestly perturbed, quasi-analytical 
model (including structures such as spiral arms or dust clumps) 
can properly recover the “true” dust mass of the input model. To 
eliminate the possibility that our fitting procedure would behave 
differently, we repeated our analysis for such basic input models, 
using exactly the same procedure as for the more realistic input 
models. We found that our fits could indeed recover the proper 
dust mass within a narrow margin. 

These results suggest that the level of dust underestimation 
is driven by the fundamental differences between the input mod¬ 
els, implying that the complex and inhomogeneous structure of 
galaxies can hide up to three times more dust than is “observed” 
when the optical images are fitted with smooth axisymmetric 
models. Although our analysis is too anecdotal to be conclusive, 
we would still argue that this effect may help explain the dust 
energy balance problem, at least in part and for certain types of 
galaxies. 
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